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Abstract. We conduct a direct comparison of three dif- 
ferent representative numerical codes for constructing 
models of rapidly rotating neutron stars in general relativ- 
ity. Our aim is to evaluate the accuracy of the codes and 
to investigate how the accuracy is affected by the choice of 
interpolation, domain of integration and equation of state. 
In all three codes, the same physical parameters, equations 
of state and interpolation method are used. We construct 
25 selected models for polytropic equations of state and 
22 models with realistic neutron star matter equations of 
state. The three codes agree well with each other (typical 
agreement is better than 0.1 % to 0.01 %) for most mod- 
els, except for the extreme assumption of uniform density 
stars. We conclude that the codes can be used for the con- 
struction of highly accurate initial data configurations for 
polytropes of index N > 0.5 (which typically correspond 
to realistic neutron stars), when the domain of integra- 
tion includes all space and for realistic equations with no 
phase transitions. With the exception of the uniform den- 
sity case, the obtained values of physical parameters for 
the models considered in this paper can be regarded as 
"standard" and we display them in detail for all models. 
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1. Introduction 

The physical state of the neutron star matter has not been 
fully understood yet because it is very difficult to inves- 
tigate particle interactions beyond nuclear matter density 
(en /c 2 ~ 2 x 10 14 g cm" 3 ) either from nuclear experiments 
or from nuclear theories, (here £n is the energy density of 
the nuclear matter and c is the velocity of light). Given 



this situation, one promising approach to explore the be- 
havior of very high density matter is to make use of the 
macroscopic quantities of neutron stars. In particular, the 
mass and the rotational period of neutron stars depend 
crucially on the softness of the equation of state (EOS) at 
very high densities (see e.g. Friedman et al. 1984, 1986, 
1989), thus, observational constrains, matched with theo- 
retical models, may help in reconstructing the equation of 
state of very high density matter. 

Given a particular equation of state, the mass of neu- 
tron stars varies with central energy density and always 
reaches a maximum. This implies that if the maximum 
mass of neutron star models constructed with a certain 
equation of state is smaller than the mass of observed 
neutron stars, that equation of state must be discarded. 
Currently, the largest accurately measured mass of slowly 
rotating neutron stars is -Mbp = 1.44M©, where Mbp is 
the mass of one of the component s of t he binary pulsar 
PSR1913+16 (Taylor & Weisberg |l989|) and M is one 
solar mass. Individual masses of neutron stars have also 
bee n est imated in si x oth er binary pulsars (Thorsett et 
al. 1993 , Wolszczan 1997), as well as in six X-ray bina- 
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ries (van Kerkwijk et al. 1995 ) but the accuracy is not as 
good as in PSR1913+16. Thus, equations of state which 
give larger masses than Mbp for slowly rotating stars, can 
be valid as candidates for the real equation of state at 
very high densities. Since the maximum mass of neutron 
stars is smaller for more compressible (soft) equations of 
state than for less compressible (stiff) equations of state, 
the true equation of state at high densities cannot be ex- 
tremely soft. 

On the other hand, stiff equations of state can be lim- 
ited by considering the neutron star with the shortest rota- 
tional period, i.e. the most rapidly rotating pulsar. There 
exists a lower limit on the rotational period for each equa- 
tion of state, because if the centrifugal force exceeds the 
self-gravity at the equatorial surface, no equilibrium states 
are allowed. The lower limit of the rotational period de- 
pends on the softness of the equation of state - the radius 



2 



T. Nozawa et al.: Highly accurate models of rotating neutron stars 



of neutron stars with softer equations of state is smaller, 
which allows for higher rotation rates. Thus, if very rapidly 
rotating neutron stars should be found, we could exclude 
most stiff equations of state. At the moment, the short- 
est period of observed pulsars is 1.56 ms, of PSR1937+21. 
Consequently, equations of state for which the shortest ro- 
tational period is larger than this value, must be excluded 
as candidates for the real equation of state for neutron 
star matter. 

The discussions above require us to make use of highly 
accurate schemes for constructing rotating neutron star 
models, in order to compute precise theoretical values of 
masses and rotational periods. Highly accurate relativistic 
equilibrium models are also needed as initial data for rel- 
ativistic time-evolution codes (modeling of nonlinear pul- 
sations, collapse and generation of gravitational waves). 
Recently, a number of groups have succeeded in construct- 
ing models of rapidly rotating neutron stars (Fried man 
et al. 1984, 1986, 1988, 1989, Eriguc hi et a l. |l994j Sal- 



gado et al. [l994a|, |l994b| , Cook et al. |l994b| , Stergioulas 



fc Friedman 1995 - for a recent review see Stergioulas 



199S ) . However, the obtained models by those authors do 



not always agree with each other (see e.g. Friedman et al. 



198S, Eriguchi et al. 1994, Salgado etal. 1994a, Cook et 



al. 1994b, Stergioulas & Friedman 1995). Although Ster- 



gioulas & Friedman (1995) have determined the cause of 
the discrepancy between models in Friedman et al. (1989) 
and Eriguchi et al. (1994), (which was due to the use of 
a slightly different equation of state table), the reasons of 
smaller differences which remain, even after using exactly 
the same equation of state, have not been clarified yet. 
This is because numerical techniques used in the different 
codes, such as the choice of parameters defining the model, 
the interpolation method, the method of integrating the 
field equations, a.s.o. are not the same. 

In this paper, three groups using their own codes (Ko- 
matsu et al. 1989a, Eriguchi et al. 1994, Salgado et al. 



1994a, Stergioulas & Friedman 1995) will decrease the dif- 



ferences between their results to a minimum possible, by 
tuning each code and using the same parameters, the same 
schemes of interpolation, the same equations of state, and 
so on. Since the basic schemes used by the three groups 
are different, it will be impossible to have exactly the same 
results and the relative differences between results are a 
measure of the accuracy of the codes. Models obtained 
with small relative differences between the three codes 
can be considered as "standard" models for each equa- 
tion of state. Furthermore, this direct comparison allows 
us to investigate the effect that the choice of interpolation 
method, equation of state and domain of integration has 
on the accuracy of the codes. 

2. Construction of neutron star models 

If rapidly rotating neutron stars were nonaxisymmetric, 
they would emit gravitational waves in a very short time 



scale and settle down to axisymmetric configurations 
Moreover, the gravity of typical neutron star is strong be- 
cause 
2GM, 



c 2 R n 



0.4, 



(1) 



where G is the gravitational constant, and M ns ~ 1.4Mq 
and R ns ~ 10km are the mass and the radius of a typical 
neutron star. Therefore, we need to solve for rotating and 
axisymmetric configurations in the framework of general 
relativity. 

For the matter and the spacetime we make the follow- 
ing assumptions: 

1. The matter distribution and the spacetime are axisym- 
metric. 

2. The matter and the spacetime are in a stationary state. 

3. The matter has no meridional motions. The only mo- 
tion of the matter is a circular one that is represented 
by the angular velocity. 

4. The angular velocity Q is constant, as seen by a distant 
observer at rest. 

5. The matter can be described as a perfect fluid. 

Under these assumptions, the metric can be expressed 
as follows (Papapetrou 1966, Carter 1969): 

ds 2 = -e 2v dt 2 + e 2a (dr 2 +r 2 d9 2 ) 

+e 2fj r 2 sm 2 d(dtp-Ludt) 2 , (2) 
= -e 2 »dt 2 +e 2 <- 2 »(dr 2 +r 2 de 2 ) 

+B 2 er 2v r 2 sin 2 9{dtp - udt) 2 , (3) 
= -e 2v dt 2 +e 2 »{dr 2 +r 2 d0 2 ) + e 2 ^{dip-ujdt) 2 , (4) 

where t is the time coordinate and the polar coordinates 
(r, 9, ip) are used. The metric depends on four metric func- 
tions (or metric potentials) which are functions of r and 
only. Different authors have used the set of functions 
(v, ui, a, (3), (u, u>, C, B)ot {v, w, /i, tp), which are related to 
each other through (g). In @ and throughout the text 
we use gravitational units (c = G = 1), unless otherwise 
stated. 

The energy momentum tensor T ab is expressed as 

T ab = (e+p)u a u b +pg ab 1 (5) 

where e, p and u a are the energy density, the pressure 
and the four- velocity, respectively. In the coordinate basis 
defined by (0), the components of the four- velocity are 



=(i,o,o,Q), 



- ( 6 ) 

where the proper velocity v with respect to a local zero 
angular momentum observer is defined by 

v = rsm9e p -"(Q-uj). (7) 

Using the metric and the energy momentum tensor 
mentioned above, we can write down the Einstein equa- 
tions for the metric components. Although we omit de- 
tailed expressions for the Einstein equations here, one can 
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easily derive them by straightforward calculations or con- 
sult the papers by Butterworth & Ipser (1976), Komatsu 
et al. (1989a) and Bonazzola et al. (1993). 

The equation of hydrostationary equilibrium can be 
derived from the equations of motion and takes the fol- 
lowing form: 



1 



1 



P 



Vp + Vv- - V ln(l - if) = 0. 



(8) 



This equation can be integrated, if we specify the equation 
of state which relates the energy density to the pressure. 
For a given EOS, a model is uniquely specified by fixing 
two parameters, such as the central energy density e c and 
the ratio of the polar to the equatorial coordinate radius, 
r p /r e . Then, four Einstein field equations and the equa- 
tion of hydrostationary equilibrium must be solved with 
appropriate boundary conditions, to yield the four met- 
ric functions and the density distribution. The available 
codes for obtaining relativistic rotating neutron star mod- 
els differ basically in the choice and method of integration 
of the four field equations and in the finite grid and finite 
difference scheme used for the integration. 



3. Codes 

We will compare three different codes, which follow 
the KEH scheme, developed by Komatsu, Eriguchi and 
Hachisu (1989a, 1989b) for general relativistic polytropes 
with uniform and differential rotation and improved by 
Cook, Shapiro & Teukolsky (1992, 1994a,b) and the 
BGSM scheme due to Bonazzola, Gourgoulhon, Salgado 
& Marck (1993). 

Concering the KEH scheme, for the present com- 
parison, we will use the original KEH code, KEH(OR), 
and the KEH code by Stergioulas & Friedman (1995), 
KEH(SF), which follows the Cook, Shapiro & Teukol- 
sky (1992, 1994a,b) modification of the KEH scheme. 

3.1. A short description of each code 

In this section the three different numerical codes are 
briefly described. Details of these codes can be found in 
Komatsu et al. (1989a, 1989b) and in Eriguchi et al. (1994) 
for the KEH(OR) code, in Stergioulas k Friedman (1995) 
for the KEH(SF) code and in Bonazzola et al. (1993) for 
the BGSM code. 



3.1.1. The KEH(OR) code 

Komatsu et al. (1989a) have developed a new scheme 
for solving rapidly rotating relativistic stars. The Ein- 
stein equations for three metric potentials v, (3 and u> are 
transformed into integral equations by using appropriate 
Green's functions for the elliptical type differential opera- 
tors. In principle, one can choose Green's functions which 



decrease as 1/r or more rapidly at large distances. Con- 
sequently boundary conditions at infinity, i.e. asymptot- 
ically flat conditions can be easily included in the inte- 
gral equations. It is noted that in this integral represen- 
tation the integrand contains the metric and the matter 
quantities such as the energy density and the pressure. 
The fourth metric a obeys a first order partial differen- 
tial equation which can be easily integrated, if the other 
metric potentials are known. The domain of integration is 
truncated at a finite distance from the star (roughly twice 
the equatorial radius) and the metric potentials in the 
integrands are assumed to vanish at that finite distance 
(instead of at infinity). 

The KEH scheme is the extended version of the self- 
consistent-field (SCF) scheme which was develope d for 



solving Newto nian rotating stars (Ostriker & Mark 1968 
Hachisu |l986p and applied to relativistic rotating stars 



by Bonazzola & Schneider (1974) with a choice of metric 
functions different from that of the codes considered in 
this article. In the SCF method, the iteration proceeds as 
follows. If one assumes initial guesses for the matter quan- 
tities and the metric potentials, new (and better) values 
for v, (3 and u> can be obtained using the integral represen- 
tations for the metric potentials The fourth metric poten- 
tial a can be easily solved as mentioned before. By using 
newly obtained metric potentials, a new density and a 
new pressure can be computed from the hydrostationary 
equilibrium equation (^). One needs to repeat the same 
procedure until the relative differences between the newly 
obtained quantities and the old ones become less than a 
certain small number, typically 10 -5 . 

In the original KEH code, the ratio of the central pres- 
sure (p c ) to the central energy density (e c ), 



Pc 



(9) 



and the ratio of the polar radius (r p ) to the equatorial ra- 
dius (r c ), r p /r , are chosen as two independent model pa- 
rameters which specify the model uniquely. The KEH(OR) 
code used in this comparison differs from the code used in 
Komatsu et al. (1989a) only in one aspect, that is men- 
tioned in the next section. 

3.1.2. The KEH(SF) code 

The KEH(SF) code (Stergioulas & Friedman, 1995) differs 
from the original KEH scheme in two ways. First, it follows 
Cook et al. (1992) in using a redefined radial variable 



(10) 



r + r e 



where r e is the value of the coordinate r at the equator. 
By this transformation, the region [0, oo] in the r coordi- 
nate is transformed to the region [0, 1] in the s coordinate. 
Consequently, the domain of integration of the field equa- 
tions does not have to be truncated at a finite distance 
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from the star (as in the KEH(OR) code) but covers all 
space. With this choice, the boundary conditions can be 
applied accurately at infinity. 

Second, Stergioulas & Friedman found that the choice 
of coordinates in the original KEH scheme results in the 
metric potential a oscillating in the radial direction. The 
oscillation is especially pronounced inside the star and in- 
troduces an error of 1-2 % in the mass, radius and other 
quantities. The problem was fixed by using a finite differ- 
ence formula for the second order radial derivative that 
uses twice the grid-spacing. Although this formula is, in 
principle, of lower accuracy, the oscillations are damped 
completely, resulting in a more accurate stellar model. 

The KEH(OR) code used in this comparison has been 
modified so as to use the same second order derivative 
formula as Stergioulas & Friedman, to smooth out the 
oscillations in the metric potential a. 

3.1.3. The BGSM code 

Bonazzola et al. (1993) have developed a new formulation 
based on the 3 + 1 formalism which has been used in hy- 
drodynamics in general relativity. Their choice of slicing 
and gauge in the 3+1 formalism results in the same form 
of the metric usually chosen for stationary problems, i.e. 

. Consequently the Einstein equations are reduced to 
the same differential equations of elliptic type as used by 
other schemes. 

The main part of the BGSM formulation is similar 
to that of the KEH formulation, except for the metric 
coefficient £ = a + v for which a second-order (elliptic) 
equation is used instead of a first-order equation in KEH. 
The Einstein equations are reorganized so as to "pick" out 
the Laplacian operators in two and three dimensional flat 
spaces and regard all other remaining terms as "source 
terms" in the Poisson-like equations in two and three di- 
mensional flat spaces. Concerning the matter, essentially 
the same equation as (||) is used for the hydrostationary 
equation. 

The characteristic features of the BGSM code can be 
found in the numerical solving method, i.e. the pseudo- 
spectral method (Gottlieb & Orszag 1977, Bonazzola et 
al. 1996, 1997, 1998b). In the spectral method all func- 
tions are expanded in terms of certain base functions and 
algebraic equations for coefficients which appear in the ex- 
pansion are solved. Therefore there are two distinct pro- 
cedures in this method: one is obtaining coefficients from 
the functions and the other is constructing functions by 
using the coefficients. Since these two steps need to be, in 
general, performed many times, it is highly desirable to 
use a fast algorithm. In the spectral method of the BGSM 
code, Bonazzola et al. (1993) have adopted trigonometric 
functions for the angle variable and the Chebyshev poly- 
nomials for the radial variable. Consequently for the angle 
part of any function, the fast fourier transform (FFT) can 
be employed. Concerning the radial variable, a new vari- 



able which is related to the radial variable by a simple 
equation is introduced so that the Chebyshev polynomi- 
als are expressed by the trigonometric functions. After this 
transformation, one can use the FFT also for the radial 
variable. 

The BGSM code can handle the region extended to 
infinity as is done by the KEH(SF) code. This can be per- 
formed by introducing a new radial variable u as follows: 

R 

u = — , for outside of the matter (11) 
r 

where Rq is the maximum value of the stellar radius. The 
boundary conditions at infinity, are easily applied at u — 
0. 

It may be fair to note that in the Newtonian rotating 
star problems a similar expansion was used by Ostriker 
and Mark (1968) in the SCF method, although Ostriker 
and Mark used the integral form of the Newtonian poten- 
tial instead of solving the Poisson equation directly and 
they did not use the FFT. 

3.2. Relations among the three different codes 

Here we summarize similarities and difference between the 
three codes: 

A) Common features through all three codes: 

1. Quasi-isotropic coordinates are used. It implies that 
the metric components are essentially the same, al- 
though background views deriving the metric are not 
the same. 

2. Integral form of the hydrostationary equation for the 
matter is used. 

3. Poisson-like operators are "picked" out from the Ein- 
stein equations and the other terms are treated as 
"source terms" . 

B) The KEH(OR) code differs from the other two 
codes in that the boundary conditions are not applied at 
infinity, but approximate boundary conditions are applied 
at a finite distance from the star. 

C) A difference between the BGSM code and the other 
codes is the use of a second-order (elliptic)equation for £ 
in BGSM versus a first-order equation for a = £ — v in 
KEH. 

We reorganized our codes so as to make the differences 
as small as possible. The codes agree exactly on: 

1. the physical model parameters by which we can specify 
a model uniquely, 

2. the values of physical constants, and 

3. the equation of state of matter. 

Since the three codes use different grids and/or nu- 
merical methods for solving the field equations, there will 
always be a residual difference in the results, even after 
this reorganization. This residual difference is what we 
want to determine. 
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3. 3. Starting point of computations 

3.3.1. Constants 

Values used in this paper for the velocity of light c, the 
gravitational constant G, the mass of the sun M & and the 
baryon mass tob are as follows: 

c = 2.9979 x 10 10 cm s _1 , 
G = 6.6732 x ltrV 1 cm 3 s" 2 , 
M Q = 1.987 x 10 33 g, 
m B = 1.66 x 10~ 24 g. 

Note that some of the above constants differ slightly 
from the ones used in the papers where the three codes 
were first presented. 

3.3.2. Interpolation of Tabulated Equation-of-State Data 

For realistic equations of state, the energy density, pres- 
sure and other thermodynamical quantities are given in ta- 
bles. Intermediate values need to be obtained by a method 
of interpolation. We will use two different interpolation 
schemes, the four-point Lagrange interpolation (hereafter 
LI) and the cubic Hermite interpolation (HI) (Swesty 
1996): 

A) Lagrange interpolation 

Let us assume that there is a table which relates the 
variable x to the variable y at n points, i.e. a set of values 
(xi,yi) for i = 1, ... ,7i are tabulated. For the LI scheme 
the interpolated formula j/li can be expressed as 



where 



p'{x)\ x = Xi ' 



p(x) = (x-xi)(x-x 2 )...(x-x n ), 
p(x) 



Pi(x) 



(x - Xi) ' 



(12) 

(13) 
(14) 



The prime denotes the differentiation with respect to the 
argument. This scheme does not guarantee the values of 
derivatives at the points in the table. In this paper we use 
the four point Lagrange interpolation, i.e. n = 4. 
B) Hermite interpolation 

In the Hermite interpolation, the interpolated formula 
j/hi for xi < x < Xi + i is expressed as 



Um{x) = yih (w) + y i+1 h (l - to) 

'dy_ 

dx 

'dy_ 

dx 



(xi+i ~ Xi)hi(w) 

(x i+1 - Xi)hi(l - w), 



i+1 



where 



to 



Xi-\-l Xi 



(15) 



(16) 



Here h Q and hi are the cubic Hermite functions defined 

by 



ho(w) = 2w 3 
hi(w) = w 3 - 



- 3w z 
2w 2 - 



1, 



to. 



(17) 
(18) 



It is noted that these cubic Hermite functions are unique 
polynomials of degree three satisfying the following rela- 
tions: 



MO) = 1, 

Mi) = Mo) = Mi) 
Mo) = i, 

Mo) = Mi) = Mi) 



o. 



(19) 
(20) 
(21) 
(22) 



Contrary to the LI, in the Hermite interpolation the values 
as well as their first derivatives at mesh points are exactly 
reproduced by the interpolation formula. The main ad- 
vantage of the Hermite interpolation is that it preserves 
the thermodynamical consistency (Swesty 1996). 

4. Equations of state 

4-1. Relativistic polytropes 

We use the fo llowin g relation as a polytropic equation of 
state (Tooper 1965): 



K- 



P' 



pc 



P 



7 



1 



(23) 
(24) 
(25) 



7-1 
Kp'< , 
1 

N ' 

where K and N are the polytropic constant and polytropic 
index, respectively, while p is the rest mass density. 

It should be noted that this equation of state includes 
the limiting case of e = pc 2 = constant, when 7 = 00 
(N = 0). The constant density models are also called 
homogeneous models. For polytropes of index N < 1.0, 
the density does not go to zero smoothly at the surface 
and the first derivatives of the density across the surface 
are discontinuous. This kind of discontinuity may become 
the cause of unfavorable behavior of solutions, unless it is 
treated carefully. For constant density models, the situa- 
tion is even worse, since the density itself is discontinuous 
across the surface. 

Although polytropic EOSs are not as realistic as tab- 
ulated EOSs (but one can reproduce neutron star bulk 
properties with an N ~ 1.0 polytrope), they are help- 
ful to check numerical codes. Since the hydrostationary 
equation can be analytically integrated and no additional 
numerical errors arise in solving it. 



4-2. Short description of realistic equations of state 

As discussed in the introduction, the main uncertain- 
ties about neutron star properties are related to the un- 
known interactions of the neutron star matter at high 
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density regions. In the last decades, many equations of 
state have been proposed by considering different kinds 
of interactions into account. A large collection of repre- 
sentative equations of state were compiled by Arnett & 
Bowers (1977), who constructed nonrotating neutron star 
models and obtained physical quantities for slowly ro- 
tating neutron stars. We will choose three equations of 
state of Arnett & Bowers' compilation, i.e. equations C, 
G and L according to their notation. Equations C, G and 
L are those derived by Bethe & Johnson (1974), Canuto 
& Chitre (1974) and Pan dhari pande & Smith (1975) (see 
also Pandharipande et al. 1976; ), respectively. Those equa- 
tions of state were also used by Friedman et al. (1986) 
for constructing rapidly rotating relativistic neutron stars 
models. 

In addition to these equations of state we also em- 
ploy the WFF3 (UV14+TNI) equation of state by Wiringa 
et al. (1988), the FPS equation of state by Lorenz et 
al. (1993), and the equation of state which represents a 
causal limit (CLES). 

Some characteristic features of each equation of state 
can be summarized as follows. 



Bethe - Johnson I (C): EOS C is of intermediate stiffness. 
The maximum gravitational mass of a spherical neutron 
star for this EOS is 1.85M Q . The density range is from 
1.71 x 10 14 g cm" 3 to 3.23 x 10 15 g cm" 3 . Hyperons as 
well as nucleons are taken into account. The interaction is 
assumed non-relativistic and represented by the modified 
Reid soft core potential with non-integer parameters. To 
include the many-body theory, the constrained variational 
principle is employed. This equation of state is joined to 
the composite BBP(e/c 2 > 4.3 x 10 n g cm" 3 ) - BPS(10 4 
g cm- 3 < e/c 2 < 4.3 x 10 n g cm" 3 ) - FMT(e/c 2 < 10 4 
g cm" 3 ). Here BBP, BPS and FMT denote equations of 
state by Baym et al. (1971a), Baym et al. (1971b) and 
Feynman et al. (1949), respectively. 

Canuto - Chitre (G): EOS G an extremely soft equation 
of state. The maximum gravitational mass of a spheri- 
cal neutron star for this equation of state is 1.36M , so 
this EOS is not acceptable as a realistic candidate for the 
true EOS of neutron star matter. It is used in this com- 
parison, because it is close to the softest possible realistic 
EOS consistent with observational constraints. The den- 
sity range is from 2.37 x 10 15 g cm" 3 to 7.23 x 10 15 g 
cm" 3 . Crystallization of neutrons is included. The inter- 
action is non-relativistic and represented by the modified 
Reid soft core potential. This equation of state is joined to 
the composite PC(7 x 10 14 g cm" 3 < e/c 2 < 2.4 x 10 15 g 
cm" 3 ) - BBP(4.3 x 10 11 g cm" 3 < e/c 2 < 7 x 10 14 g 
cm" 3 ) - BPS(10 4 g cm" 3 < e/c 2 < 4.3 x 10 n g cm" 3 ) - 
FMT(e/c 2 < 10 4 g cm -3 ). Here PC denotes the equation 
of state by Pandharipande (1971). 



Pandharipande - Smith (L): EOS L is an extremely stiff 
EOS. The maximum gravitational mass of a spherical neu- 
tron star for this equation of state is 2.70M Q . The den- 
sity range is larger than 4.386 x 10 11 g cm" 3 . Composi- 
tions consist of neutrons. The interaction is assumed non- 
relativistic and is represented by the nuclear attraction 
due to scalar particle exchange. This equation of state is 
joined to the BPS (10 4 g cm -3 < e/c 2 < 4.4 x 10 n g cm" 3 ) 
- FMT(e/c 2 < 10 4 g cm" 3 ). 

Wiringa - Fiks - Farbrocini (WFF3): EOS WFF3 
(Wiringa et al., 1988) is of intermediate stiffness. At 
present, the WFF3 equation of state is regarded as one of 
the best candidates for the high density region. This EOS 
is an improved version of the equation of state by Fried- 
man & Pandharipande (1981). The nucleon-nucleon inter- 
action described by a two-body Urbana UV14 potential 
and the phenomenological three-nucleon TNI interaction 
are taken into account. Compositions are considered to be 
neutrons. The maximum gravitational mass of a spherical 
neutron star for this equation of state is 1.84M . Although 
the usual WFF3 EOS is joined to EOS NV (Negele & Vau- 
thcrin, 1973), we will also use a different version, in which 
it is joined to the more modern FPS EOS (Lorenz et al., 
1993). 

Lorenz - Ravenhall - Pethick (FPS): This EOS is also 
a modern version of the equation of state by Friedman 
& Pandharipande (1981). The nucleon-nucleon interaction 
described by a two-body Urbana UV14 potential and the 
phenomenological three-nucleon TNI interaction are taken 
into account. In the FPS equation of state the Skyrme 
model is used, where the effective interaction has the spa- 
tial character of a two-body delta function plus deriva- 
tives. The FPS equation of state can be considered to be 
an improved version of the BBP equation of state in the 
region of the lower density. 

Causal limit equation of state (CLES): As an extreme 
case, we consider an equation of state which consists of a 
causal limit EOS (e = p + constant) for e/c 2 > 1.66 x 10 14 
g cm -3 and the FPS EOS below that density. The causal 
limit EOS has the property that, in the interior of the 
star, the phase velocity of sound is equal to the velocity 
of light in vacuo, i.e. v s = ^/dp/de = c. 

5. Results 

5.1. Models selected for poly tropes 

We have started our comparison project by selecting sev- 
eral representative polytropic models, the parameters of 
which are shown in Tables |l|-|[ We have chosen 

1. models of very low central density (nearly Newtonian) 
with slow and rapid rotation, 
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Table 1. Polytropic models 



Model 


N05sn 


N05rn 


N05sr 


N05mr 


N05rr 


N075sn 


N075rn 


N075sr 


N075mr 


N 


0.5 


0.5 


0.5 


0.5 


0.5 


0.75 


0.75 


0.75 


0.75 


£c 


1.00e-03 


1.00e-03 


1.20e+00 


1.04e+00 


2.00c+00 


3.00e-05 


3.00e-05 


8.00e-01 


6.50e-01 


r p /r e 


9.80e-01 


5.99e-01 


9.85e-01 


5.50e-01 


7.57e-01 


9.80e-01 


6.82e-01 


9.79e-01 


5.66e-01 


D, 


6.91e-03 


2.85e-02 


2.14e-01 


9.78e-01 


1.00e+00 


1.10e-03 


4.11e-03 


1.77e-01 


6.30e-01 




6.92e-03 




2.19e-01 


9.80e-01 








1.79e-01 


6.31e-01 


M 


6.30e-08 


9.43e-08 


1.26e-01 


1.53e-01 


1.35e-01 


3.73e-07 


4.84e-07 


1.40e-01 


1.64e-01 








1.27e-01 


1.54e-01 


1.36e-01 






1.41e-01 


1.65e-01 


Mo 


6.30e-08 


9.43e-08 


1.52e-01 


1.83e-01 


1.59e-01 


3.73e-07 


4.84e-07 


1.60e-01 


1.87e-01 








1.53e-01 




1.60e-01 










Rcirc 


3.04e-02 


4.20e-02 


4.03e-01 


5.35e-01 


4.03e-01 


1.95e-01 


2.47e-01 


5.33e-01 


7.37e-01 




















7.38e-01 


J 


1.31e-13 


1.47e-12 


2.16e-03 


1.72e-02 


1.07e-02 


4.55e-12 


3.35e-ll 


2.87e-03 


1.75e-02 




1.32e-13 




2.21e-03 






4.57e-12 




2.89e-03 




I 


1.90e-ll 


5.17e-ll 


1.01e-02 


1.76e-02 


1.07e-02 


4.16e-09 


8.15e-09 


1.62e-02 


2.77e-02 


T/W 


5.18e-03 


1.26e-01 


5.01e-03 


1.47e-01 


8.76e-02 


4.93e-03 


8.96e-02 


5.94e-03 


l.lle-01 




5.20e-03 




5.15e-03 


1.49e-01 


8.96e-02 


4.96e-03 




5.95e-03 


1.12e-01 


Zp 


2.10e-06 


3.10e-06 


6.32e-01 


9.64e-01 


9.97e-01 


1.94e-06 


2.55e-06 


4.59e-01 


6.10e-01 








6.42e-01 


9.75e-01 


10.0e-01 




2.56e-06 


4.61e-01 


6.13e-01 




-2.08e-04 


-1.19e-03 


4.10e-01 


-3.62e-01 


-1.42e-01 


-2.12e-04 


-1.01e-03 


2.61e-01 


-3.19e-01 


-2.09e-04 




4.15e-01 


-3.64e-01 


-1.47e-01 


-2.13e-04 




2.65e-01 


-3.20e-01 




2.12e-04 


1.20e-03 


8.66e-01 


2.99e+00 


2.86e+00 


2.16e-04 


1.02e-03 


6.59e-01 


1.72e+00 




2.13e-04 




8.90e-01 


3.02e+00 


2.92e+00 


2.17e-04 




6.70e-01 


1.73e+00 


e 


1.99e-01 


8.01e-01 


1.31e-01 


6.95e-01 


5.64e-01 


1.96e-01 


7.32e-01 


1.66e-01 


7.02e-01 




2.05e-01 




1.72e-01 


7.02e-01 


5.74e-01 


2.09e-01 




1.97e-01 


7.06e-01 



2. models of high central density (relativistic) with slow 
and rapid rotation, and 

3. models at the maximum mass for each EOS. 

Note that, the maximum mass model almost coincides 
with the maximum angular velocity model, unless there 
is a large phase transition at densities close to the cen- 
tral density of the maximum mass star (cf. Cook, Shapiro 
& Teukolsky, 1994b and Stergioulas & Friedman, 1995). 
For all EOSs in this comparison the two models almost 
coincide. 

In order to evaluate the performance of our numerical 
codes for models with discontinuous density distribution, 
we also compare a number of homogeneous models which 
cover both highly relativistic and Newtonian, rapidly ro- 
tating and nonrotating cases, as shown in Tables || and [| 
(the contents of these tables will be described in Sect. 5.4). 



5.3.1. Grid and physical quantities 



5.2. Models selected for realistic equations of state 

As discussed in the previous section, we use six represen- 
tative realistic EOSs: C, G, L, WFF3+FPS, WFF3+NV 
and FPS. In addition, we use the causal limit EOS CLES. 
For each equation of state, we compute several models as 
shown Tables ^| _ @- The models correspond to the max- 
imum mass model, a fast rotating 1.4A/ Q model and a 
nonrotating model for each EOS. 



For this comparison project, KEH(OR) and KEH(SF) 
have used grids with (71 x 201) and (261 x 401) 
(angular x radial) grid-points. In the equatorial plane, half 
of the radial grid-points are inside the star. BGSM uses 21 
x 41 or 33 x 65 grid points (note that the notion of "grid 
points" is not very significant for a spectral method; the 
above numbers should better be referred to as the num- 
bers of basis functions employed in the expansions of the 
physical fields). 



5.3. Computed Quantities 



Here we summarize the notation of computed physical 
quantities: 
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Table 2. Polytropic models (continued) 



Model 


NIOsn 


NIOrn 


NIOsr 


NIOmr 


NIOrr 


N15sn 


N15rn 


N15sr 


N15mr 


N15rr 


N 


1.0 


1.0 


1.0 


1.0 


1.0 


1.5 


1.5 


1.5 


1.5 


1.5 


£c 


1.00e-06 


1.00e-06 


4.00e-01 


3.40e-01 


1.00e+00 


1.00e-09 


1.00e-09 


6.50e-01 


6.10e-02 


1.50e-01 


r p /r e 


9.76e-01 


6.39e-01 


9.72e-01 


5.84e-01 


8.34e-01 


9.60e-01 


7.08e-01 


8.68e-01 


6.20e-01 


8.40e-01 


Q 


2.00e-04 


7.00e-04 


1.26e-01 


3.77e-01 


4.00e-01 


6.53e-06 


1.58e-05 


1.61e-01 


l.lle-01 


1.20e-01 












4.01e-01 












M 


2.54e-06 


3.29e-06 


1.65e-01 


1.88e-01 


1.61e-01 


9.76e-05 


l.lle-04 


2.10e-01 


2.91e-01 


2.66e-01 








1.66e-01 


1.89e-01 














M 


2.54e-06 


3.29e-06 


1.82e-01 


2.07e-01 


1.73e-01 


9.76e-05 


l.lle-04 


2.08e-01 


3.04e-01 


2.77e-01 


Rcirc 


1.27e+00 


1.71e+00 


7.82e-01 


1.09e+00 


6.79e-01 


5.29e+01 


6.63e+01 


1.27e+00 


2.85e+00 


1.79e+00 


J 


2.14e-10 


1.50e-09 


4.24e-03 


2.02e-02 


9.48e-03 


3.60e-07 


1.26e-06 


1.02e-02 


3.88e-02 


2.13e-02 








4.25e-03 
















I 


1.07e-06 


2.14e-06 


3.36e-02 


5.36e-02 


2.37e-02 


5.51e-02 


7.99e-02 


6.32e-02 


3.50e-01 


1.78e-01 










5.37e-02 






8.00e-02 


6.33e-02 






T/W 


5.56e-03 


9.10e-02 


6.61e-03 


8.36e-02 


3.48e-02 


7.48e-03 


5.24e-02 


1.32e-02 


4.75e-02 


2.29e-02 




5.57e-03 




6.63e-03 


8.41e-02 


3.52e-02 








4.76e-02 


2.30e-02 


Z p 


2.04e-06 


2.71e-06 


3.23e-01 


4.04e-01 


4.56e-01 


1.91e-06 


2.25e-06 


2.57e-01 


1.94e-01 


2.29e-01 




2.05e-06 


2.74e-06 


3.26e-01 


4.05e-01 


4.58e-01 




2.27e-06 


2.58e-01 


1.95e-01 


2.30e-01 




-2.52e-04 


-1.19e-03 


1.54e-01 


-2.83e-01 


-5.95e-02 


-3.44e-04 


-1.05e-03 


-5.33e-02 


-2.24e-01 


-8.10e-02 






1.56e-01 


-2.84e-01 


-6.19e-02 






-5.43e-02 


-2.25e-01 


-8.17e-02 


7 b 


2.56e-04 


1.20e-03 


4.97e-01 


1.15e+00 


1.03e+00 


3.48e-04 


1.05e-03 


5.75e-01 


6.23e-01 


5.49e-01 


2.57e-04 




5.02e-01 




1.04e+00 






5.78e-01 




5.51e-01 


e 


2.16e-01 


7.70e-01 


2.03e-01 


7.14e-01 


4.75e-01 


2.80e-01 


7.06e-01 


4.44e-01 


7.30e-01 


4.93e-01 




2.27e-01 




2.25e-01 


7.18e-01 


4.85e-01 


2.84e-01 


7.07e-01 


4.52e-01 




5.00e-01 



e c Central energy density 

r p /r e Ratio of polar to equatorial radii 

ft Angular velocity of the star 

Mo Baryon mass 

M p Proper mass 

M Gravitational mass 

Rdrc Equatorial circumferential radius 

r e Equatorial coordinate radius 

J Total angular momentum 

/ Moment of inertia about the rotation axis 

T Rotational energy 

W Gravitational energy 

v e Velocity of comoving observer at the equator 

relative to the locally non-rotating observer 

Z p Polar redshift 

Z c Central redshift 

Z\ q Equatorial redshift in the backward direction 

Z[ q Equatorial redshift in the forward direction 

e Intrinsic eccentricity of the star's surface 

GRV2 Two dimensional virial identity 

GRV3 Three dimensional virial identity 

Some of the quantities in the above list can be ex- 
pressed as follows: 

Z p = e- 2v ? - 1, 



7 b _ 



1-Vc 
1+V ( 

l + v t 
l-v r 



1/2 



1 + r e e(^-A)/2, 



- 1, 



1/2 



1 - r e e^-MI\ 



(26) 
(27) 

(28) 



where subscripts p and e denote values at the pole and the 
equatorial surface, respectively. 



M = 2?r / p 



M v = 2tt / e 



M = 2tt 



e 2a+/3 

. r 2 sinQdrdO, 

s 2q+/3 

vr^2 



r sin 8drd8, 



(29) 
(30) 



2p 



1 — ir 



1-V 2 

r 2 sin 9drd9 



J = 2tt / e 2 "+ 2 ^ V :"T r 3 sin 2 0drd0, 



1-v 2 



(31) 
(32) 



T = - I VidJ = 2tt / e 



W = Mpc 2 + T - Mc 2 
and 



2a+20 (g+P)« o. .; ; , 2 
1-V 2 



fir 6 sin^ 0drd6, (33) 
(34) 

(35) 



The eccentricity of the meridional cross section is dc 
fined by the following procedure (Friedman et al. 1986 ) 
If the surface of the star is defined by 



r = r s (9), 



(36) 
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the metric of the stellar surface can be expressed as 



Table 3. Spherical constant density models 



da 2 = e 2 V sin 2 8dip 2 + e 



2<i 



dr s 



r s (8? 



de 2 . (37) 



If we embed this surface in the flat three dimensional 
space, it is expressed as 



R = R s (z), 



(38) 



in cylindrical coordinates (R, <p, z). The 2-metric of this 
surface is 



da 2 



dR s 

dz 



+ 1 



dz 2 + R 2 s dip 2 



(39) 



Comparing these two equations, we have the following re- 
lations, if they express the same surface geometry: 



Rs(B) 

and 

Zs(0) -- 



e'V sin#, 








r • \ 


e 2a 



(40) 



dR s 

de 



1/2 



Using these quantities the eccentricity is defined as 



e=\ 1 



z s (e = o) 

Rs{e = 7T/2) 



(41) 



(42) 



For polytropes, it is convenient to express quantities 
in dimensionless form, by using 

K N/2 

as a fundamental 

length scale, as was done in Cook et al. (1994a). In ge- 
ometrized units (c = G = 1), dimensionless quantities are 
define as follows: 



R 



r = 




(43) 


circ — 


R- N/2 R ■ 


(44) 




K N ^ 2 n, 


(45) 


e ee 


K N e, 


(46) 


P = 


K N p, 


(47) 


P 


K N p, 


(48) 


./ EE 


K~ N J, 


(49) 


/ EE 


K- 3N / 2 J, 


(50) 


M EE 


K- N ' 2 M, 


(51) 


M EE 


K- N ' 2 M a . 


(52) 



5.3.2. Virial Theorem 

Equilibrium configurations in Newtonian gravity satisfy 
the following relation: 



Model 


NOOsn 


NOOsr 


£c 


1.00e+00 


1.00e+00 


Pc 


1.00e-04 


1.00e+00 


r p /r e 


1.00e+00 


1.00e+00 


M 


1.38e-06 


1.12e-01 




1.41e-06 


1.13e-01 


Mo 


1.38e-06 


1.59e-01 




1.41e-06 


1.61e-01 


Rcirc 


6.91e-03 


2.99e-01 




6.93e-03 


3.00e-01 




2.00e-04 


9.71e-01 




2.04e-04 


10.1e-01 


7f 
^eq 


2.00e-04 


9.71e-01 


2.04e-04 


10.1e-01 


rrb 


2.00e-04 


9.71e-01 




2.04e-04 


10.1e-01 


Table 4. 


Rotating constant de 


Model 


NOOrn 


NOOrr 


£c 


1.00e+00 


1.00e+00 


Pc 


1.00e-04 


1.00e+00 


r p /r e 


6.50e-01 


7.00e-01 


n 


1.00e+00 


1.40e+00 




1.02e+00 


1.41e+00 


M 


2.05e-06 


1.35e-01 




2.06e-06 


1.39e-01 


M 


2.05e-06 


1.84e-01 




2.06e-06 


1.87e-01 


Rcirc 


9.06e-03 


3.45e-01 




9.11e-03 


3.46e-01 


J 


6.76e-ll 


1.37e-02 




6.97e-ll 


1.41e-02 


I 


6.74e-ll 


9.82e-03 




6.87e-ll 


10.0e-03 


T/W 


1.08e-01 


1.63e-01 




1.14e-01 


1.68e-01 


Zp 


2.83e-04 


1.60e+00 




2.86e-04 


1.71e+00 


Zlq 


-8.81e-03 


-1.55e-01 


-8.96e-03 


-1.60e-01 


7 b 

^eq 


9.38e-03 


9.41e+00 


9.54e-03 


11.3e+00 


e 


7.60e-01 


7.07e-01 




7.61e-01 


7.11e-01 



where U is the internal energy. This relation is called the 
virial relation and has been used to check the accuracy of 
numerically obtained solutions (see e.g. Ostriker & Mark 



1968, Tassoul 1978) 



2T + 3( 7 - l)U + W = 0, 



(53) 



In general relativity, similar relations were first found 
by Bonazzola (1973). Recently, two virial identities in gen- 
eral relativity have been discovered by Gourgoulhon & 
Bonazzola (1994) and Bonazzola & Gourgoulhon (1994). 
Those identities are valid for a general asymptotically flat 
spacetime. We can use these identities to estimate the nu- 
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Table 5. Spherical models with realistic EOSs 



Model 


Csr 


Gsr 


Lsr 


WFF(FPS)sr 


WFF(NV)sr 


FPSsr 


CLESsr 


EOS 


C 


G 


L 


WFF3+FPS 


WFF3+NV 


FPS 


CLES 


£ c [g cm^] 


1.09e+15 


6.31e+15 


4.30e+14 


1.22e+15 


1.22e+15 


1.31e+15 


1.85e+14 


r p /r e 


1.00e+00 


1.00e+00 


1.00e+00 


1.00e+00 


1.00e+00 


1.00e+00 


1.00e+00 


M[M Q ] 


1.41E+00 


1.36E+00 


1.39e+00 


1.41e+00 


1.41e+00 


1.41e+00 


1.41e+00 






1.37E+00 


1.41e+00 










M [M Q ] 


1.55e+00 


1.57c+00 


1.49c+00 


1.57e+00 


1.57c+00 


1.57c+00 


1.51c+00 








1.53e+00 










Rdrc [km] 


1.19e+01 


6.94e+00 


1.48e+01 


1.09e+01 


1.09e+01 


1.08e+01 


1.77e+01 


Zp 


2.39e-01 


5.28e-01 


1.76e-01 


2.69e-01 


2.68e-01 


2.72e-01 


1.43e-01 




2.41e-01 


5.38e-01 


1.80e-01 


2.71e-01 


2.71e-01 


2.74e-01 


1.44e-01 




2.39e-01 


5.28e-01 


1.76e-01 


2.69e-01 


2.68e-01 


2.72e-01 


1.43e-01 


2.41e-01 


5.38e-01 


1.80e-01 


2.71e-01 


2.71e-01 


2.74e-01 


1.44e-01 


7 b 


2.39e-01 


5.28e-01 


1.76e-01 


2.69e-01 


2.68e-01 


2.72e-01 


1.43e-01 


2.41e-01 


5.38e-01 


1.80e-01 


2.71e-01 


2.71e-01 


2.74e-01 


1.44e-01 



Table 6. Rotating models with realistic EOSs 



Model 


Cbr 


Cmr 


Gbr 


Gmr 


Lbr 


L(L)mr 


L(H)mr 


WFF(FPS)br 


EOS 


C 


C 


G 


G 


L 


L(LI) 


L(HI) 


WFF3+FPS 


e c [g cm - ' 5 ] 


8.70e+14 


2.64e+15 


3.10e+15 


5.58e+15 


3.90e+14 


1.20e+15 


1.20c+15 


9.70e+14 


r p /r e 


6.73e-01 


5.72e-01 


6.37e-01 


5.77e-01 


7.08e-01 


5.53e-01 


5.53e-01 


6.30e-01 


n[io 4 s- 1 ] 


5.89e-01 


1.07e+00 


1.16c+00 


1.57e+00 


4.21e-01 


8.17e-01 


8.13e-01 


7.00e-01 










1.58e+00 


4.25e-01 


8.18e-01 


8.14e-01 


7.01e-01 


M[M e ) 


1.41e+00 


2.17e+00 


1.41e+00 


1.57e+00 


1.41e+00 


3.31e+00 


3.29e+00 


1.42e+00 












1.43e+00 


3.33e+00 


3.30e+00 




M [M Q ] 


1.52e+00 


2.48e+00 


1.58e+00 


1.79e+00 


1.50e+00 


3.90e+00 


3.87e+00 


1.55e+00 




1.53e+00 


2.49e+00 


1.59c+00 




1.53e+00 


3.92e+00 


3.88e+00 




_R circ [km] 


1.52c+01 


1.32e+01 


1.02c+01 


9.24e+00 


1.76c+01 


1.85e+01 


1.85e+01 


1.43e+01 






1.34e+01 




9.25e+00 


1.77e+01 








cJ/GMl 


1.20e+00 


2.99e+00 


1.19e+00 


1.54e+00 


1.24e+00 


7.72e+00 


7.62c+00 


1.31e+00 






3.00e+00 






1.28c+00 


7.80e+00 


7.64c+00 




/[10 45 g cm 2 ] 


1.80e+00 


2.45e+00 


9.01e-01 


8.56e-01 


2.60e+00 


8.31e+00 


8.23e+00 


1.64e+00 






2.46c+00 


9.03e-01 


8.59e-01 


2.64c+00 


8.40e+00 


8.25c+00 


1.65e+00 


T/W 


8.44e-02 


1.10e-01 


9.59e-02 


1.05e-01 


8.38e-02 


1.37e-01 


1.36e-01 


1.01e-01 




8.49e-02 




9.69e-02 


1.07e-01 


8.47c-02 


1.39e-01 




1.02e-01 


Zp 


2.46e-01 


6.86e-01 


4.63e-01 


7.44e-01 


1.90e-01 


8.26e-01 


8.18e-01 


2.82e-01 






6.89e-01 


4.66e-01 


7.49e-01 


1.93e-01 


8.31e-01 


8.21e-01 


2.84e-01 


^eq 


-1.81e-01 


-3.30e-01 


-2.42e-01 


-3.36e-01 


-1.42e-01 


-3.47e-01 


-3.46e-01 


-2.10e-01 




-3.31e-01 


-2.43e-01 


-3.38e-01 


-1.43e-01 


-3.49e-01 


-3.47e-01 


-2.11e-01 


^eq 


6.96e-01 


1.94c+00 


1.28e+00 


2.10e+00 


5.44e-01 


2.43e+00 


2.41c+00 


8.11e-01 




6.98e-01 


1.96c+00 




2.12e+00 




2.46e+00 


2.42e+00 


8.12e-01 


e 


6.84e-01 


6.90e-01 


6.75e-01 


6.74e-01 


6.65e-01 


6.97e-01 


7.02c-01 


7.13e-01 




6.88e-01 


7.26e-01 


6.82e-01 


6.81e-01 


6.70e-01 


7.04e-01 


7.03e-01 


7.18e-01 



merical error. Let us define two quantities A2 and A3 as 
follows: 



A 2 = 



and 



-OO />7T 



JO 



p+(e+p) 



1-v 2 



e 2fl rdrd9 



Jo 



dvdv-'-e^-^dudu 
4 



-. (54) 



rdr d9 



e 2 ^rdrd9 




dvdv — -d/j,dtp 



1 d/j. 



\dr r tan 9 89 , 

dip 1 dtp 



+ l(l-e 2 ^r 2 sm 2 9)(^ + 
4r V or 



1 



r sin 2 9) 8 



dudtu 



rt&n9 89 
\ -1 



e^rdrd9 
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Table 7. Rotating models with realistic EOSs (continued) 



Model 


WFF(FPS)mr 


WFF(NV)br 


WFF(NV)mr 


FPSbr 


FPSmr 


CLESbr 


CLESmr 








(max. mass) 




(max. mass) 




(max. mass) 


EOS 


WFF3+FPS 


WFF3+NV 


WFF3+NV 


FPS 


FPS 


CLES 


CLES 


£ c [g cm - ' 5 ] 


2.70e+15 


9.70e+14 


2.71e+15 


1.02e+15 


2.91e+15 


1.81e+14 


4.20e+14 


r p /r e 


5.65e-01 


6.30e-01 


5.65e-01 


6.40e-01 


5.68e-01 


5.90e-01 


5.31e-01 


fi[10 4 s- 1 ] 


1.15e+00 


6.98e-01 


1.15e+00 


6.97e-01 


1.18c+00 


3.63e-01 


6.03e-01 






6.99e-01 










6.06e-01 


M[M Q ] 


2.19e+00 


1.41e+00 


2.19e+00 


1.41e+00 


2.12c+00 


1.41e+00 


6.64e+00 




2.20c+00 


1.42c+00 


2.20e+00 




2.13e+00 


1.42e+00 


6.71e+00 


M O [M ] 


2.55e+00 


1.54e+00 


2.55e+00 


1.54e+00 


2.46e+00 


1.50e+00 


8.36e+00 






1.55e+00 










8.40e+00 


R circ [km] 


1.28c+01 


1.43e+01 


1.28e+01 


1.41e+01 


1.24e+01 


2.30e+01 


2.86e+01 
















2.87e+01 


cJ/GMl 


3.22e+00 


1.31e+00 


3.21e+00 


1.27c+00 


2.96e+00 


1.62e+00 


3.61e+01 




3.23c+00 




3.22e+00 




2.97e+00 




3.64e+01 


/[10 45 g cm 2 ] 


2.45e+00 


1.64e+00 


2.44c+00 


1.60c+00 


2.20e+00 


3.92e+00 


5.25e+01 




2.46c+00 


1.65c+00 


2.46e+00 




2.21e+00 


3.93e+00 


5.30e+01 


T/W 


1.22e-01 


1.00e-01 


1.22e-01 


9.70e-02 


1.17e-01 


1.13e-01 


1.92e-01 




1.23e-01 


1.01e-01 


1.23e-01 


9.76e-02 


1.18e-01 




1.97e-01 


Zp 


7.63e-01 


2.81c-01 


7.62e-01 


2.81e-01 


7.47e-01 


1.60e-01 


1.46c+00 




7.64e-01 


2.83e-01 


7.64e-01 


2.82e-01 


7.55e-01 


1.61e-01 


1.49e+00 


Zeq 


-3.39e-01 


-2.10e-01 


-3.38e-01 


-2.05e-01 


-3.37e-01 


-1.92e-01 


-4.02e-01 


-3.41e-01 


-2.11e-01 


-3.40e-01 


-2.06e-01 


-3.39e-01 


-1.93e-01 


-4.07e-01 




2.19e+00 


8.10e-01 


2.18e+00 


8.03e-01 


2.15c+00 


5.24e-01 


5.89e+00 


2.21c+00 


8.12e-01 


2.20e+00 


8.05e-01 


2.17e+00 


5.25e-01 


6.14e+00 


e 


6.93e-01 


7.13e-01 


6.93e-01 


7.05e-01 


6.84e-01 


7.69e-01 


7.20e-01 




7.28e-01 


7.18e-01 


7.27e-01 


7.10e-01 


6.90e-01 


7.70e-01 


7.29e-01 



with the abridged notation 
dfidip 1 d/j, dtp 

Then, we define: 

GRV2 = |1-A 2 |, 
GRV3 = |1-A 3 |. 



(56) 

(57) 
(58) 



If the Einstein equations are satisfied, these quantities sat- 
isfy the following virial identities: 



GRV2 = 0, 
GRV3 = 0. 



(59) 
(60) 



Since exact solutions for the stationary problems satisfy 
the above relations, we can choose GRV2 and GRV3 as the 
error indicators for numerically obtained solutions. Note 
that due to its three dimensional character, GRV3 gives a 
larger weight to the external layers of the star. GRV3 is a 
relativistic generalization of the Newtonian virial theorem 
Eq. (43). 

In practice, however, it should be noted that the virial 
identities in the above form are not always close to the 
accuracy of numerical results. In particular, for GRV2 the 
integration is done by integrating in the 9 coordinate as 
seen from the definition of GRV2. If one does not use the 
9 coordinate for solving equilibrium structures, one needs 
to change variables and in that procedure accuracy may 



be lost and the resultant values may become worse than 
the "real" accuracy before the variable change. 

Concerning the quantity GRV3, the metric potentials 
in the vacuum region contribute to the integral consider- 
ably. It implies that if only the finite regions are treated, 
as in the KEH(OR) code, a large portion of the integrand 
cannot be taken into account, However, the expressions 
for the GRV2 and GRV3 are not unique because we can 
do the integral by part and replace the second derivatives 
with the matter terms, using Einstein's equations. In this 
way, the contribution far away from the star becomes less 
important. In the KEH(OR) code, GRV2 and GRV3 are 
evaluated through 




p+(e + p) 



(1-v 2 ) 



e 2a rdrd9 



V, r + —V, { 



1 1 

-v \ u, rr +-v, r +^v,ee 



rdr d9 



3 ,r max ,7T 

4 Jo Jo 



Id 



167re 



2a (fi-w)(e+p) 



l-v 2 



1 COt 9 ,ns 
"-W, r -I —UJ,0 +UJ, r [p + U), r 



12 



T. Nozawa et al.: Highly accurate models of rotating neutron stars 



~U>,0 (P + V),0 



(61) 



3p+ (e + p) 




e 2a+f3 r 2 



, ^ + O0ap)- 2rtanedr 

+ HE (i _ e 2 "- 2 / 3 ) 

4r tan 9 dr 



(1"« 2 ) 
-v (d 2 v + dvdj3) 
1 da 



1-e 



2a-2f3\ 



-—Lul6ire 



rV sin 9 drd9 



2:i 



_ 2v {e+ P )(n-u) 
(l-v 2 ) 



e 2 ^ 2 V sin 2 i 



with another abridged notation 



(62) 



0V 



a 2 v> 2 ^ i a 2 v> 



<9r 2 



<9r 



+ 



1 



dtfj 



86 1 



r 2 tan 6 89 



(63) 



and r max is the distance of the truncated point beyond 
which actual numerical computations are not carried out 
in the KEH(OR) code. 

This rewrite does not break the mathematical iden- 
tity. In a sense, it may make "identity" more trivial, and 
then what information the identities provide us becomes 
unclear. However, as far as the same expression of the 
identity in the same code is used, they can play a role as 
indicators of accuracy among models solved by each code. 

5.4- Tables of Models and Comparison 

The physical parameters of 45 models, computed in this 
comparison project, are displayed in Tables 0^0- All 
quantities are displayed to three significant figures. The 
lower and the upper bounds on each quantity, as obtained 
by comparing results from the three codes, are shown in 
the upper and lower rows for each corresponding quantity, 
respectively. It follows that quantities expressed in a sin- 
gle row can be regarded as "exact", to three significant 
figures. 

Tables [j] and || display results for polytropes with in- 
dex N — 0.5, 0.75, 1.0 and 1.5. For each value of the 
polytropic index N we compute the following models: 

1. a spherical Newtonian model (denoted by the symbol 

sn), 

2. a rapidly rotating Newtonian model (rn), 

3. a nearly spherical, relativistic model (sr) , 

4. the maximum mass model (mr) and 

5. a rapidly rotating relativistic model (rr). 



For the constant density case (N = 0), the spheri- 
cal Newtonian and spherical relativistic models are dis- 
played in Table || and the rapidly rotating Newtonian and 
rapidly rotating relativistic models in Table [l| While all 
other models are specified by the central energy density 
e c and the ratio of the polar radius to the equatorial ra- 
dius r p /r e , constant density models are specified by their 
central pressure and r p /r e . 

For realistic equations of state, spherical models are 
shown in Table || and rotating models with 1.4M Q (br: 
binary pulsar mass and relativistic) and maximum mass 
models are shown in Tables || and [?]. For the equation of 
state L, model L(L)mr uses four-point Lagrange interpo- 
lation, while L(H)mr is the same model but computed 
using cubic Hermite interpolation. 

From the tables displaying polytropic models one can 
see that the three codes have a good agreement on most 
quantities especially for soft polytropes. For stiff poly- 
tropes (N < 1.0) the agreement is somewhat smaller. For 
constant density models, the relative differences between 
the three codes become several percent. More sensitive 
quantities are the three redshifts and the eccentricity. It 
should be noted that redshift factors are local quantities 
which reflect the metric potentials at each point. This im- 
plies that local values of the metric potentials do not have 
the same agreement between different numerical codes, 
as integrated global quantities. For the eccentricity, one 
needs to compute the length along the surface of the star 
(see the definition of the eccentricity ((42|)). Since in the 
KEH codes fi = cos 9 is used as the angular variable there 
arise numerical errors near the pole region, i.e. 9 ~ 0. 
Thus, the differences in the values of the eccentricity also 
reflect this numerical error due to the choice of coordi- 
nates. This causes differences of up to a few percent in 
the eccentricity for rapidly rotating models. On the other 
hand, global quantities such as angular velocity, mass, ra- 
dius and angular momentum agree quite well among re- 
sults of different codes. 

From Tables || - [?], similar tendencies can be observed 
for realistic equations of state. Models for the most EOSs, 
except EOS CLES, have a good agreement between the 
three codes, although the agreement is not as good as 
for polytropic models.. By comparing models constructed 
with EOSs WFF(FPS) and WFF(NV), it is evident that 
the choice of the low density EOSs affects very little the 
structure of the star. 

The main reason for the large differences in the con- 
stant density case is that the discontinuous density dis- 
tribution is creating Gibbs phenomena near the surface 
and this affects all three codes. The reason for the smaller 
agreement for realistic EOSs, compared to polytropic 
EOSs, is that the necessary interpolation between tab- 
ulated data affects the accuracy with which the equation 
of hydrostationary equilibrium is satisfied. For EOS L, the 
choice of the interpolation scheme also affects the accuracy 
of the computed models, with the cubic Hermite scheme 
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being a better choice compared to a four-point Lagrange 
interpolation (see the discussion in a later section) . For the 
other realistic EOSs the choice of the interpolation scheme 
had a negligible effect on the accuracy of computed mod- 
els. 



5.5. Detailed comparison 

In order to investigate further the differences among nu- 
merical results obtained by the three codes, we show more 
detailed results for models: N15sn, N15mr and N05sn in 
Table |, N05mr in Table |, NOOsn, NOOrn and NOOsr 
in Table NOOrr in Table ||, Gmr, Lsr and L(L)mr 
in Table ||, L(H)mr, WFF(FPS)sr and WFF(FPS)6r in 
Table [l3| and CLESmr in Table In these tables, values 
to eight figures for each physical quantity are shown, as 
well as the relative differences among results of the three 
codes. 

The three relative differences diffl, diff2 and diff3 are 
defined as 



diffl 
diff2 
diff3 



BGSM — KEH(OR) 
= KEH(OR) ' 

KEH(OR) — KEH(SF) 
= KEH(SF) 

KEH(SF) - BGSM 
= BGSM ' 



(64) 
(65) 
(66) 



From Tables || and ^, we see that the agreement be- 
tween the KEH(SF) and BGSM codes for the rapidly ro- 
tating, relativistic polytropic models models N15mr and 
N05mr is between 10 -4 and 3x 10~ 4 in all computed quan- 
tities (global and local) , except for the eccentricity (due to 
reasons we explained in the previous section). This very 
good agreement shows that both the BGSM and KEH 
schemes (the latter when applying boundary conditions 
exactly at infinity) are suitable for the construction of 
highly accurate initial-data configurations of rapidly ro- 
tating relativistic stars, modeled as polytropes with typi- 
cal index N > 0.5. 

From the same tables, we see that the agreement be- 
tween the KEH(OR) and BGSM codes is similar to the 
agreement between KEH(SF) and BGSM in the global 
quantities of model N15mr but to within 10 -3 for the 
local quantities of this model. For the stiffer polytrope 
N05mr the agreement is 10~ 3 and 10~ 2 for global and 
local quantities, respectively. This difference in accuracy 
between KEH(SF) and KEH(OR) is expected, since in 
KEH(OR) boundary conditions are applied only approxi- 
mately at a finite distance from the star. Considering how 
close to the star the domain of integration is truncated, 
the KEH(OR) code performs very well. This is explained 
as follows: 



Since the integration is performed over only a finite 
region, the truncated part of the integral, I tr , can be ex- 
pressed as 



S{v )G(r,r )d 3 r 



(67) 



where S and G are the source term and the Green's func- 
tion, respectively, and the subscript out denotes that the 
integration covers the region with r > R m ax- Here R ma x 
is the maximum radius of the region of integration. Al- 
though the source terms of the integrals for the metric 
potentials in the KEH(OR) scheme contain both matter 
and the metric terms, only the metric terms contribute 
to the integral the outside of the star. Since we impose 
asymptotically flat conditions for the metric potentials, 
the radial dependence of the source term can be consid- 
ered to be 



S 



1 



(68) 



where k > 4, because of the behavior of metric functions 
at large radii. Consequently, I tr can be roughly expressed 
as 



Jt r (r) ~ constant, 
1 



In 1 - 



[T <C r e <C Rmax 

r 



Rr, 



(r e <r< R max ) (69) 



From the above equation, we can see the following: First, 
for the stellar interior, i.e. r < r e , Itr ~ constant, so that 
the relative error It r /Io can be very small where Jo is the 
exact value of the integral. Second, for the outer region 
of the stars, i.e. r e < r < R m ax, the value of Itr becomes 
larger for larger values of r because of the logarithmic 
term. Third, for the same region, i.e. r e < r < R m ax, if the 
source term depends weakly on the radial coordinate, i.e., 
for smaller values of k, the contribution of the truncated 
part becomes larger. This occurs for the potentials v and 
oj. However, for the metric functions B and £, since values 
of k are large, the relative differences are rather small. 

In the extreme case of constant density models (Ta- 
bles [lO] and |ll|) , the three codes agree on the computed 
physical quantities typically only within a few percent and 
this is caused by the sharp density discontinuity at the sur- 
face of the star. The numerical schemes in this comparison 
assume that the density distribution is a smooth function 
of coordinates, thus, in the case of density discontinuities, 
this assumption is violated and Gibbs phenomena appear, 
resulting in low accuracy of the computed models. 

From Tables p^!4] it follows that the agreement of the 
KEH(SF) code to the BGSM code is between 1CT 3 and 
1CT 4 for realistic EOSs, except for EOS CLES, where the 
agreement is an order of magnitude smaller. KEH(OR) 
and BGSM agree on the realistic models within 10~ 2 and 
10 -3 , i.e. similar to the agreement for the N — 0.5 poly- 
trope. The somewhat lesser agreement for realistic EOSs 
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Table 8. Detailed comparison of polytropic models 



Model 


KEH(OR) 


KEH(SF) 


BGSM 


diffl 


difI2 


difT3 


Model 


N15sn 












e c 


1.0000000e-09 


1.0000000e-09 


1.0000000e-09 








r p /r e 


9.6000000e-01 


9.6000000e-01 


9.6000000e-01 








0, 


6.5347464e-06 


6.5349925e-06 


6.5346091e-06 


-2e-5 


-4e-5 


6e-5 


M 


9.7631159e-05 


9.7633565e-05 


9.7634578e-05 


4e-5 


-3e-5 


-le-5 


Mo 


9.7631237e-05 


9.7633644e-05 


9.7634657e-05 


4e-5 


-3e-5 


-le-5 


Rcirc 


5.2907355e+01 


5.2907099e+01 


5.2908476e+01 


2e-5 


5e-6 


-3e-5 


J 


3.5985713e-07 


3.5988285e-07 


3.5988120e-07 


7e-5 


-7e-5 


5e-6 


I 


5.5068262e-02 


5.5070123e-02 


5.5073103e-02 


9e-5 


-3e-5 


-5e-5 


T/W 


7.4801796e-03 


7.4805723e-03 


7.4802183e-03 


5e-6 


-5e-5 


5e-5 




1.9104721e-06 


1.9105286e-06 


1.9104977e-06 


le-5 


-3e-5 


2e-5 




-3.4382698e-04 


-3.4383282e-04 


-3.4382701e-04 


9e-8 


-2e-5 


2e-5 




3.4764792e-04 


3.4765387e-04 


3.4764801e-4 


3e-7 


-2e-5 


2e-5 


e 


2.7990046e-01 


2.8386885e-01 


2.8014870e-01 


9e-4 


-le-2 


le-2 


GRV2 


1.0578017e-04 


6.3653021e-05 


-9.0173678e-08 








GRV3 


-1.3712946e-05 


7.8924860e-05 


-2.2726769e-07 








Model 


N15mr 












e c 


6.1000000e-02 


6.1000000e-02 


6.1000000e-02 








r p /r e 


6.2000000e-01 


6.2000000e-01 


6.2000000e-01 








Q 


1.1082917e-01 


1.1080248e-01 


1.1079616e-01 


-3e-4 


2e-4 


6e-5 


M 


2.9099968e-01 


2.9091548e-01 


2.9091876e-01 


-3e-4 


3e-4 


-le-5 


M 


3.0433314e-01 


3.0436572e-01 


3.0437123e-01 


le-4 


-le-4 


-2e-5 


Rcirc 


2.8538213e+00 


2.8537335e+00 


2.8539175e+00 


3e-5 


3e-5 


-6e-5 


J 


3.8794961e-02 


3.8804941e-02 


3.8806669e-02 


3e-4 


-3e-4 


-5e-5 


I 


3.5004286e-01 


3.5021727e-01 


3.5025283e-01 


6e-4 


-5e-4 


-le-4 


T/W 


4.7633988e-02 


4.7508685e-02 


4.7507826e-02 


-3e-3 


3e-3 


2e-5 


Zp 


1.9429844e-01 


1.9478524e-01 


1.9478308e-01 


3e-3 


-3e-3 


le-5 


Z[q 


-2.2475260e-01 


-2.2447647e-01 


-2.2448160e-01 


-le-3 


le-3 


-2e-5 


7 b 
^eq 


6.2263043e-01 


6.2336170e-01 


6.2335930e-01 


le-3 


-le-3 


4e-6 


e 


7.2977397e-01 


7.3040218e-01 


7.3099356e-01 


2e-3 


-9e-4 


-8e-4 


GRV2 


3.8566454e-04 


3.8935202e-03 


-4.5069715e-07 








GRV3 


-5.3198333e-05 


1.1087078e-04 


-2.5949035e-06 








Model 


N05sn 












£c 


1.0000000e-03 


1.0000000e-03 


1.0000000e-03 








r p /r e 


9.8000000e-01 


9.8000000e-01 


9.8000003e-01 








a 


6.9205090e-03 


6.9211147e-03 


6.9078733e-03 


-2e-3 


-9e-5 


2e-3 


M 


6.2992939e-08 


6.2995192e-08 


6.2989922e-08 


-5e-5 


-4e-5 


8e-5 


M 


6.2993012e-08 


6.2995264e-08 


6.2989995e-08 


-5e-5 


-4e-5 


8e-5 


Rcirc 


3.0433836e-02 


3.0433775e-02 


3.0432960e-02 


-3e-5 


2.e-6 


3e-5 


J 


1.3150934e-13 


1.3153761e-13 


1.3125388e-13 


-2e-3 


-2e-4 


2e-3 


I 


1.9002842e-ll 


1.9005263e-ll 


1.9000620e-ll 


-le-4 


-le-4 


2e-4 


T/W 


5.1994036e-03 


5.2005153e-03 


5.1792787e-03 


-3e-3 


-2e-4 


4e-3 


Zp 


2.0985378e-06 


2.0986158e-06 


2.0984061e-06 


-6e-5 


-3e-5 


le-4 


Zeq 


-2.0851997e-04 


-2.0851108e-04 


-2.0812950e-04 


-2e-3 


4e-5 


2e-3 


7 b 
^eq 


2.1271705e-04 


2.1270832e-04 


2.1232632e-04 


-2e-3 


4e-5 


2e-3 


e 


1.9884365e-01 


2.0484010e-01 


1.9911770e-01 


le-3 


-3e-2 


3e-2 


GRV2 


-7.8558660e-06 


2.5493102e-05 


-8.1720158e-06 








GRV3 


-1.1614383e-05 


9.1053508e-05 


-2.6069983e-06 









is due to the use of interpolation between the tabulated 
equation of state data (see the discussion next). 

In Tables || to [l4|, we also display the virial quantities 
GRV2 and GRV3. In the ideal case, these should exactly 
vanish, so the smaller the values for GRV2 and GRV3 
are, the better is the accuracy of the computed model. 
The opposite is not always true, i.e. in some models the 



computed values for GRV2 or GRV3 do not reflect an 
overall better agreement in physical quantities among the 
different codes. This indicates that the computation of 
GRV2 and GRV3 may itself be prone to numerical error. 
This seems to be the case for GRV2 in rapidly rotating 
models computed with the KEH(SF) code, where one first 
has to interpolate data between cos 9 and 9 grids to be able 
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Table 9. Detailed comparison of polytropic models (continued) 



Model 


KEH(OR) 


KEH(SF) 


BGSM 


diffl 


diff2 


diff3 


Model 


N05mr 












e c 


1.0400000e+00 


1.0400000e+00 


1.0400000e+00 








r p /r e 


5.5000000e-01 


5.5000000e-01 


5.5000156e-01 








ft 


9.8004159e-01 


9.7793068e-01 


9.7787248e-01 


-2e-3 


2e-3 


6e-5 


AI 


1.5376721e-01 


1.5302143e-01 


1.5301513e-01 


-5e-3 


5e-3 


4e-5 


M 


1.8250033e-01 


1.8251809e-01 


1.8250776e-01 


4e-5 


-le-4 


6e-5 


Rcirc 


5.3491240e-01 


5.3492344e-01 


5.3489664e-01 


-3e-5 


-2e-5 


5e-5 


J 


1.7214828e-02 


1.7214474e-02 


1.7217664e-02 


2e-4 


2e-5 


-2e-4 


I 


1.7565405e-02 


1.7602959e-02 


1.7607269e-02 


2e-3 


-2e-3 


-2e-4 


T/W 


1.4949334e-01 


1.4721775e-01 


1.4726559e-01 


-2e-2 


2e-2 


-3e-4 


Zp 


9.6404034e-01 


9.7533189e-01 


9.7519945e-01 


le-2 


-le-2 


le-4 




-3.6445066e-01 


-3.6225810e-01 


-3.6208301e-01 


-7e-3 


6e-3 


4e-4 


yb 


2.9851556e+00 


3.0239590e+00 


3.0243455e+00 


le-2 


-le-2 


-le-4 


e 


6.9783608e-01 


6.9495486e-01 


7.0174527e-01 


6e-3 


4e-3 


-le-2 


GRV2 


-2.9470467e-03 


3.5351133e-02 


-5.3131831e-05 








GRV3 


-7.2928314e-02 


1.1781519e-04 


1.3585788e-04 









to compute GRV2. Note that the displayed values of the 
two virial quantities for the KEHYOR) code correspond 
to the modified virial identities (|6l]) and (|6^) and not to 
the original identities. The computation of GRV3 for the 
KEH(OR) code is affected significantly by the truncation 
of the domain of integration. 

6. Discussion and Conclusion 

6.1. Discussion 

6.1.1. Metric Potentials 

As redshift factors differ by about 10 % for the constant 
density, relativistic model NOOrr between the three codes 
(while the agreement global quantities is within a few %) 
we compare directly the local values of metric potentials 
for several models. Figures 1 to 16 show the four met- 
ric potentials (upper panel) and the relative differences in 
them between BGSM and KEH(OR) (middle panel) and 
KEH(SF) and BGSM (lower panel) for the models N05mr, 
N15mr, L(L)mr and WFF(FPS)6r. The metric potentials 
are graphed against the coordinate r in the equatorial 
plane (8 — tt/2, solid line) and along the axis of rota- 
tion (9 = 0, dashed line). The range of the coordinate r is 
the twice the equatorial radius of the star. 

In general, the agreement in the local values of the 
metric potentials reflects the agreement in the computed 
physical parameters of models. In these graphs, several 
significant behaviors can be pointed out: First, there are 
high frequency and small amplitude oscillations at the in- 
ner part of the stars for all models. Second, the differences 
between the results of KEH(OR) and those of the other 
two codes are growing outside the stars as r increases. 
Third, although the differences between the KEH(SF) and 
BGSM codes are very small for models N05mr, N15mr 
and WFF(FPS)6r, there appear larger differences for the 
stiff model L(L)mr. Fourth, there appears a larger ampli- 



tude oscillation in the metric potential lj on the axis of 
rotation, close to the surface. 

The first behavior is due to the integration scheme 
of the KEH code, i.e. the Simpson scheme. In general, 
the Simpson scheme gives results with higher precision, 
compared with those obtained by the trapezoidal scheme. 
However, in the KEH scheme, the integrands contain non- 
smooth functions with respect to the radial coordinate, 
because of the nature of the Green's functions. The non- 
uniform distribution of the weight factor in Simpson's 
scheme for nonsmooth functions results in oscillating be- 
haviors with very small amplitudes, which cannot be no- 
ticed in the behavior the integrated quantities. 

The second behavior in the original KEH code is 
caused by the "truncation" of the domain of integration 
at a finite distance from the star, instead of integrating 
over the whole space. 

The large differences in the metric potentials between 
KEH(SF) and BGSM for EOS L, could be accounted to 
the stiffness of the equation of state, but the differences 
between KEH(OR) and BGSM for the same model are not 
as large, and we have not an explanation for that. 

The oscillations in oj on the axis of rotation near the 
surface are present also for the soft N = 1.5 polytropes, 
while for N = 0.5 they are larger. This indicates that 
terms in the field equations for lo are very sensitive to 
the presence of the surface and the accompanying Gibbs 
phenomenon. Even for N = 1.5 polytropes, where the 
density goes to zero smoothly at the surface, there is a 
small scale Gibbs phenomenon, due to the finite number 
of grid points used to represent the region of integration. 

6.1.2. Method of Interpolation 

An important factor for the local accuracy of models con- 
structed with realistic equations of state is the method of 
interpolation of the energy vs. pressure data given in an 
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Table 10. Detailed comparison of constant density models 



Model 


KEH(OR) 


KEH(SF) 


BGSM 


diffl 


diff2 


diff3 


Model 


NOOsn 












£c 


l.OOOOOOOc+00 


1.0000000e+00 


l.OOOOOOOc+OO 








pc 


1.0000000e-04 


1.0000000e-04 


1.0000000e-04 








r p /r e 


l.OOOOOOOc+00 


l.OOOOOOOc+00 


l.OOOOOOOc+00 








M 


1.4049747e-06 


1.4079464e-06 


1.3811478e-06 


-2e-2 


-2e-3 


2e-2 


M 


1.4051472e-06 


1.4081229e-06 


1.3813135e-06 


-2e-2 


-2e-3 


2e-2 




6.9250048e-03 


6.9085570e-03 


6.9085014e-03 


-2e-3 


2e-3 


8e-6 


Zp 


2.0197813e-04 


2.0386139e-04 


1.9998000e-04 


-le-2 


-9e-3 


2e-2 




2.0197813e-04 


2.0386139e-04 


1.9998000e-04 


-le-2 


-9e-3 


2e-2 


7 b 


2.0197813e-04 


2.0386139e-04 


1.9998000e-04 


-le-2 


-9e-3 


2e-2 


GRV2 


8.5024396e-03 


2.5005801e-02 


1.9992785e-10 








GRV3 


1.1669659e-02 


3.1034530e-02 


2.4994651e-10 








Model 


NOOrn 












£c 


l.OOOOOOOc+OO 


l.OOOOOOOc+OO 


l.OOOOOOOc+00 








Pc 


1.0000000e-04 


1.0000000e-04 


1.0000000e-04 








r p /r e 


6.5000000e-01 


6.5000000e-01 


6.4996883e-01 


-5e-5 




5e-5 


Q 


1.0038565e+00 


1.0145815c+00 


1.0154202c+00 


le-2 


-le-2 


-8e-4 


M 


1.8558447e+00 


1.8466846c+00 


1.8532465e+00 


-le-3 


5e-3 


-4e-3 


M 


2.0497584e-06 


2.0637739e-06 


2.0569322e-06 


4e-3 


-7e-3 


3e-3 


M 


2.0500389e-06 


2.0640573e-06 


2.0572047e-06 


4e-3 


-7e-3 


3e-3 


Ftcirc 


9.0583317e-03 


9.1109677e-03 


9.1096124e-03 


6e-3 


-6e-3 


2e-4 


J 


6.7616039e-ll 


6.9668544e-ll 


6.9165592e-ll 


2e-2 


-3e-2 


7e-3 


I 


6.7356280e-ll 


6.8667275e-ll 


6.8115244e-ll 


le-2 


-2e-2 


8e-3 


T/W 


1.0792359e-01 


1.1 0890 17e-01 


1.1416842e-01 


6e-2 


-3e-2 


-3e-2 


Zp 


2.8331355e-04 


2.8580524c-04 


2.8505153e-04 


6e-3 


-9e-3 


3e-3 


%lq 


-8.8147112e-03 


-8.9620804e-03 


-8.9698802e-03 


2e-2 


-2e-2 


-9e-4 


7 b 
^eq 


9.3813731e-03 


9.5337271e-03 


9.5400451e-03 


2e-2 


-2e-2 


-7e-4 


e 


7.5985797e-01 


7.6083298e-01 


7.5926925e-01 


-8e-4 


-le-3 


2e-3 


GRV2 


-1.0486413e-03 


4.0913020e-02 


1.5297023e-03 








GRV3 


-1.9636423e-03 


6.3495351e-04 


1.8732802e-03 








Model 


NOOsr 












£c 


l.OOOOOOOc+OO 


l.OOOOOOOc+OO 


l.OOOOOOOc+00 








pc 


l.OOOOOOOc+OO 


l.OOOOOOOc+00 


l.OOOOOOOc+00 








r p /r e 


1.0000000e+00 


1.0000000e+00 


l.OOOOOOOc+00 








M 


1.1404379e-01 


1.1257086e-01 


1.1220252e-01 


-2e-2 


le-2 


3e-3 


Mo 


1.5912882e-01 


1.6139149e-01 


1.5914795e-01 


le-4 


-le-2 


le-2 


Rcirc 


3.0003488e-01 


2.9923839e-01 


2.9920671e-01 


-3e-3 


3e-3 


le-4 


Zp 


9.7107121e-01 


1.0138916c+00 


l.OOOOOOOc+00 


3e-2 


-4e-2 


le-2 


Z[q 


9.7107121e-01 


1.0138916c+00 


l.OOOOOOOc+00 


3e-2 


-4e-2 


le-2 


7 b 

^eq 


9.7107121e-01 


1.0138916c+00 


l.OOOOOOOc+00 


3e-2 


-4e-2 


le-2 


GRV2 


-4.2458406e-03 


8.0302649e-03 


1.2600343e-10 








GRV3 


3.4088425e-02 


7.3549407e-03 


1.3790569c-10 









EOS table. While global quantities are not affected signifi- 
cantly, the virial identities for realistic EOSs, are sensitive 
to the interpolation scheme This can be considered to re- 
flect the nature of the interpolation scheme as mentioned 
before. If we define the enthalpy (H) by 



H = ln 



e +p 
pc 2 



(70) 



the Gibbs-Duhem relation, which follows directly from the 
first law of thermodynamics, implies 



dp 

dH =£+P - 



(71) 



In the cubic Hermitc interpolation, the Gibbs-Duhem rela- 
tion is used to replace by VH the term Vp/(e + p) which 
appears in the hydrostationary equilibrium equation. If 
the tabulated function p(H) fails to satisfy the above 
relation, then the hydrostationary equilibrium equation, 
which is derived from the Bianchi identity, is only approx- 
imately verified by the numerical solution, which results 
in increased error in the GRV2 and GRV3 virial identities. 

The four point Lagrange interpolation does not satisfy 
the Gibbs-Duhem relation because it only reproduces the 
values of the discrete points, but there is no guarantee for 
the reproduction of the derivatives. This explains why the 
GRV2 and GRV3 errors are bad, even in the nonrotating 
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Table 11. Detailed comparison of constant density models (continued) 



Model 


KEH(OR) 


KEH(SF) 


BGSM 


diffl 


diff2 


diff3 


Model 


NOOrr 












£c 


1.0000000e+00 


1.0000000e+00 


1.0000000e+00 








Pc 


1.0000000e+00 


1.0000000e+00 


l.OOOOOOOc+00 








r p /r e 


7.0000000e-01 


7.0000000e-01 


7.0075459e-01 


le-3 




-le-3 


Q, 


1.3961457e+00 


1.4071531e+00 


1.3980528e+00 


le-3 


-8e-3 


7e-3 


M 


1.8323492e+00 


1.8121685e+00 


1.8158693e+00 


-9e-3 


le-2 


-2e-3 


M 


1.3889820e-01 


1.3605365e-01 


1.3462398e-01 


-3e-2 


2e-2 


le-2 


Mo 


1.8665283e-01 


1.8693714e-01 


1.8377945e-01 


-2e-2 


-2e-3 


2e-2 


Rcirc 


3.4513863e-01 


3.4566026e-01 


3.4455898e-01 


-2e-3 


-2e-3 


3e-3 


J 


1.3838952e-02 


1.4064876e-02 


1.3733153e-02 


-8e-3 


-2e-2 


2e-2 


I 


9.9122546e-03 


9.9952703e-03 


9.8230574e-03 


-9e-3 


-8e-3 


2e-2 


T/W 


1.6825846e-01 


1.6281418e-01 


1.6338670e-01 


-3e-2 


3e-2 


-4e-3 


Zp 


1.6030994e+00 


1.7071394e+00 


1.6720522e+00 


4e-2 


-6e-2 


2e-2 


z[ q 


-1.5970481e-01 


-1.5951673e-01 


-1.5541165e-01 


-3e-2 


le-3 


3e-2 


7 b 


9.4121600e+00 


1.1342393e+01 


1.0436972e+01 


le-1 


-2e-l 


9e-2 


e 


7.0676794e-01 


7.1111349e-01 


7.0812060e-01 


2e-3 


-6e-3 


4e-3 


GRV2 


-1.5336628e-02 


4.5755715e-02 


1.7034437e-03 








GRV3 


-1.0930542e-01 


9.5032627e-04 


4.2345737e-03 










o> 0.00002 
o 

1 l-e-05 

'-3 o.o 

■a -l.e-05 

" -0.00002 

-0.00003 
0.00003 

<u 0.00002 















Hiii 

" is! 1 lii 


4 








BGSM-KEH(OR) 











- K / 
_ i V 


KEH(SF)-BGSM 



H l.e-05 
'S 0.0 
B -l.e-05 

" -0.00002 
-0.00003 



Fig. 1. Metric potential B as a function of coordinate ra- 
dius for model N15mr (upper panel). Relative difference 
of B for the same model constructed with the BGSM and 
KEH(OR) codes (middle panel) and with the KEH(SF) 
and BGSM codes (lower panel). The solid graph corre- 
sponds to = 7r/2 (equatorial plane) and the dashed line 
to 9 = (axis of rotation). The largest value of r dis- 
played, corresponds to twice the coordinate radius of the 
surface of the star in the equatorial plane. 



case (GRV2 = 3E-03, GRV3 = 1E-02 for model Lsr) as 
compared to GRV2 ~ 10~ 14 for polytropic models (see 
e.g. Bonazzola et al. 



1993) 



The GRV2 and GRV3 error 
indicators thus do not reflect the precision of the code but 
the bad thermodynamical behavior of the tabulated EOS. 

The advantage of the cubic Hermite interpolation is 
that the Gibbs-Duhem relation is automatically satisfied 
because this interpolation reproduces not only the values 
themselves but also the derivatives (Swesty 1996). More- 
over, in our case, the energy density and the baryon num- 
ber density are obtained by 



e = 



p dlogp 



H dlogH 

e + p 



P- 



exp(-ff) 



(72) 
(73) 



Because of these equations, the Gibbs-Duhem relation is 
satisfied in every point. Note also that the value of e ob- 
tained in this way coincides exactly with e% at the points 
in the tabulated equation of state. 

6.2. Conclusion 

The comparison of three different codes for constructing 
rapidly rotating relativistic neutron star models demon- 
strates that the BGSM and KEH schemes used are highly 
accurate for typical polytropic models - when the field 
equations are solved to infinity - and for models con- 
structed with realistic equations of state, that do not have 
phase transitions. If one approximates neutron stars as 
constant density stars, then Gibbs phenomena at the dis- 
continuous surface reduce the accuracy of the computed 
models. If high accuracy in such models and in models 
with phase transitions is desired, then modified numeri- 
cal schemes - free of Gibbs phenomena - need to be used. 
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Table 12. Detailed comparison of models with realistic EOSs 



Model 


KEH(OR) 


KEH(SF) 


BGSM 


diffl 


diff2 


diff3 


Model 


Gmr 












£ c [g cm' 3 ] 


5.5828200e+15 


5.5828200e+15 


5.5828192e+15 








r p /r e 


5.7654000e-01 


5.7654000e-01 


5.7654388e-01 








Q [10 4 s" 1 ] 


1.5778207c+00 


1.5729910e+00 


1.5731590c+00 


-3e-3 


3e-3 


-le-4 


M [Mq] 


1.5705842e+00 


1.5670837e+00 


1.5652708e+00 


-3e-3 


2e-3 


le-3 


Mo [Mq] 


1.7915598c+00 


1.7930714e+00 


1.7913915c+00 


-9e-5 


-8e-4 


9e-4 


Rcirc [Aj77l] 

cJ/GMq 
I [10 45 5 cm 2 ] 


9.2371483c+00 


9.2454157e+00 


9.2471834c+00 


le-3 


-9e-4 


-2e-4 


1.5360080c+00 


1.5378972e+00 


1.5369015e+00 


6e-4 


-le-3 


7e-4 


8.5553062e-01 


8.5923863e-01 


8.5859068e-01 


4e-3 


-4e-3 


8e-4 


T/W 


1.0661337e-01 


1.0541507e-01 


1.0550611e-01 


-le-2 


le-2 


-9e-4 


Zp 


7.4862286e-01 


7.4468655e-01 


7.4429886e-01 


-6e-3 


5e-3 


5e-4 




-3.3808636e-01 


-3.3618145e-01 


-3.3645247e-01 


-5e-3 


6e-3 


-8e-4 


7 b 


2.1006841c+00 


2.1151171e+00 


2.1140605c+00 


6e-3 


-7e-3 


5e-4 


e 


6.7374530e-01 


6.8063333e-01 


6.7954077e-01 


9e-3 


-le-2 


2e-3 


GRV2 


1.3388534e-03 


1.7877189e-02 


-3.6339736e-04 








GRV3 


-3.2386612e-02 


9.7104987e-04 


-5.2521686e-04 








Model 


Lsr 












e c [g cm,- 3 ] 


4.2995000e+14 


4.2995000e+14 


4.2995391e+14 








r p /r e 


l.OOOOOOOc+00 


l.OOOOOOOc+00 


l.OOOOOOOc+00 








M [Mq] 


1.4097969e+00 


1.3884591e+00 


1.4100000c+00 


le-4 


2e-2 


-2e-2 


M [Mq] 


1.5235811c+00 


1.4938211e+00 


1.5253919c+00 


le-3 


2e-2 


-2e-2 


Rcirc \k77l] 


1.4784119e+01 


1.4819571e+01 


1.4778437c+01 


-4e-4 


-2e-3 


3e-3 


Zp 


1.7822863e-01 


1.7557838e-01 


1.7974175e-01 


9e-3 


2e-2 


-2e-2 


Zeq 


1.7829440e-01 


1.7557838e-01 


1.7974175e-01 


8e-3 


2e-2 


-2e-2 


7 b 

^eq 


1.7829440e-01 


1.7557838e-01 


1.7974175e-01 


8e-3 


2e-2 


-2e-2 


GRV2 


-3.9884976e-03 


1.7391300e-02 


3.1884440e-03 








GRV3 


2.4064437e-03 


2.6180956e-02 


1.2340895e-02 








Model 


L(L)mr 












e c [g cm' 3 ] 


1.2022600e+15 


1.2022600e+15 


1.2022594e+15 








r p /r e 


5.5300000e-01 


5.5300000e-01 


5.5300318e-01 








n [10 4 s- 1 ] 


8.1797409e-01 


8.1681214e-01 


8.1674028e-01 


-2e-3 


le-3 


9e-5 


M [Mq] 


3.3254178c+00 


3.3067743e+00 


3.3156417c+00 


-3e-3 


6e-3 


-3e-3 


Mo [Mq] 


3.9164459e+00 


3.8979008e+00 


3.9209347e+00 


le-3 


5e-3 


-6e-3 


cJ/GMq 


1.8472119e+01 


1.8460489e+01 


1.8475482e+01 


2e-4 


6e-4 


-8e-4 


7.7807280e+00 


7.7201999e+00 


7.8026429e+00 


3e-3 


8e-3 


-le-2 


I [10 45 g cm 2 ] 


8.3595026e+00 


8.3065196e+00 


8.3959624e+00 


4e-3 


6e-3 


-le-2 


T/W 


1.3918072e-01 


1.3669628e-01 
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2e-2 
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Zp 
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-3.4699020e-01 
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4e-3 
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-6e-3 


-5e-3 


e 
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7.0355618e-01 


7.0273768e-01 


8e-3 


-9e-3 


le-3 


GRV2 


9.4730085e-04 


3.4941975e-02 


-2.2424315e-03 








GRV3 


-4.0789432e-02 


2.4558001e-04 


-2.1163857e-03 









Such numerical schemes could employ, for example, sur- 
face fitted coordinates. Such a scheme has been presented 
recently by Bonazzola et al. (1998a) in the framework of 
spectral methods and looks promising for rotating stellar 
models. Further, we demonstrated that the metric poten- 
tials are subject to various local oscillatory behaviors, even 
if integrated quantities have very good accuracy. This ob- 
servation is important for the effort of constructing initial 
data for the numerical evolution of rotating relativistic 
neutron star models. 
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Fig. 2. Same as Figure 1 but for the metric potential B 
of model WFF(FPS)6r. 
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Fig. 4. Same as Figure 1 but for the metric potential B 
of model L(L)mr. 
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Fig. 3. Same as Figure 1 but for the metric potential B 
of model N05mr. 



Fig. 5. Same as Figure 1 but for the metric potential v of 
model N15mr. 
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Fig. 6. Same as Figure 1 but for the metric potential v of 
model WFF(FPS)&r. 



Fig. 8. Same as Figure 1 but for the metric potential v of 
model L(L)mr. 
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Fig. 7. Same as Figure 1 but for the metric potential v of 
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Fig. 9. Same as Figure 1 but for the metric potential u> Fig. 10. Same as Figure 1 but for the metric potential u> 
(in units of Q) of model N15mr. (in units of f2) of model WFF(FPS)6r. 
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Fig. 11. Same as Figure 1 but for the metric potential u> 
(in units of J7) of model N05mr. 



Fig. 12. Same as Figure 1 but for the metric potential u> 
(in units of f2) of model L(L)mr. 
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Fig. 13. Same as Figure 1 but for the metric potential ( 
of model N15mr. 
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Fig. 14. Same as Figure 1 but for the metric potential £ Fig. 16. Same as Figure 1 but for the metric potential £ 
of model WFF(FPS)6r. of model L(L)mr. 
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